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LARGE  AMPLITUDE  WAVES  GENERATED 


BY  AN  EXPLOSIVE  MECHANISM 


1 .  INTRODUCTION 


During  the  period  1946  to  1958  a  series  of  observations  were  made  of 
waves  formed  by  detonating  nuclear  devices  below  the  ocean  surface.  Because 
the  effects  of  deep  underwater  nuclear  explosion  are  largely  of  military 
interest,  phenomena  associated  with  these  tests  are  not  well  documented  in  the 
open  literature. 


For  a  nuclear  device  detonated  in  deep  water,  the  wave  intensity  is 
represented  by  the  peak -to -peak  height  H  of  the  wave  envelope.  The  manner  in 
which  the  height  diminishes  with  distance  or  radius  R  from  surface  zero  has 
empirically  been  found  to  be  (Glasstone  and  Dolan  (1977)): 

H~  4.36  x  10V5V\  (1.1) 


where  H  and  R  are  in  metres  and  the  yield  w  in  kilotons  of  TNT.  The  constant 
of  proportionality  in  (1.1)  diminishes  as  the  water  depth  decreases  and 
consequently  the  energy  delivered  to  the  water  is  greater  at  depth  than  for  a 
shallow  burst  having  the  same  yield.  The  empirical  expression  holds  provided 
the  depth  of  water  dw  (metres),  in  which  the  surface  waves  are  produced  is  in 
the  range: 
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(1.2) 


where  the  lower  limit  is  the  maximum  diameter  of  the  gas  bubble.  The 
relation  is  valid  for  any  depth  of  water  in  which  an  explosive  device  is 
detonated.  The  idealized  surface  displacement  is  shown  in  figure  1. 
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FIGURE  1  An  explosion  generated  wave  train  as 
observed  at  a  given  distance  from ' 
surface  zero  (Glasstone  and  Dolan  (1977)). 


A  theoretical  interpretation  of  the  above  wave  structure  presents 
some  difficulties  and  this  arises  from  the  fact  that  to  give  a  complete 
description  of  the  phenomena  would  require  details  of  the  detonation, 
propagation  of  the  shock  and  pulsating  bubble  as  well  as  venting  and  blast 
effects.  To  include  all  these  physical  phenomena  into  a  single  set  of 
equations  in  order  to  delineate  the  associated  wave  train  would  render  the 
problem  intractable  analytically. 

In  this  paper  a  mathematical  model  is  developed  which  encompasses 
(i)  the  observations  of  waves  generated  at  the  Bikini  nuclear  test  site  in  the 
Pacific  Ocean  (Glasstone  and  Dolan  (1977)),  and  (ii)  the  wave  structure  that 
can  be  expected  from  a  submarine  volcanic  explosion.  We  assume  that  the  waves 
generated  by  the  above  phenonema  are  similar  and  can  be  treated  by  the  same 
model.  Moreover  we  consider  a  rigid  cyclindrical  volume  of  water  to  be 
projected  vertically  from  a  depth  h  ,  for  infinite  water  above  a  rigid  bottom 
with  a  known  velocity  under  the  action  of  ordinary  forces.  The  method  of 
approach  will  be  to  derive  analytically  a  solution  of  the  linearized  wave 
equations  for  water  of  uniform  depth,  in  the  context  of  the  linear  theory  of 
long  waves  (Whitham  (1974)). 

The  above  type  of  wave  disturbance  can  also  be  attributed  to 
tectonic  plate  uplifting  of  oceanic  blocks  leading  to  the  formation  of  tsumani 
waves.  Such  problems  are  of  interest  in  dynamical  oceanography. 

The  waves  generated  by  the  above  phenomena  are  also  of  interest  in 
understanding  the  vulnerability  of  ocean  as  well  as  coastal  structures  to  wave 
loading  and  inundation. 


EQUATIONS  OF  MOTION 


j 
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Consider  an  inviscid  incompressible  fluid  in  a  constant 
gravitational  field.  A  rigid  cylinder  of  fluid  of  radius  0  <  r  <  rQ  is 
assumed  to  be  projected  upwards  from  a  depth  z  «  -hD,  to  the  plane  of  the 
undisturbed  free  surface  z  *  0  with  a  velocity  £(t)  in  the  positive  z 
direction,  as  shown  in  figure  2.  It  will  be  assumed  that  each  element  of 
moving  fluid  experiences  no  nett  rotation  from  one  instant  to  the  next. 
Therefore  the  fluid  is  assumed  to  be  irrotational,  i.e.  there  is  zero 
vorticity  in  the  fluid.  If  u  is  the  fluid  velocity  then  curlu  *  0.  For 
irrotational  flow  we  may  write 

u  =  V*,  (2.1 ) 

where  $  is  the  velocity  potential  and  7  the  gradient  operator.  For 
incompressible  flow: 

V.Ur.«  0,-.  (2.2) 

and  substituting  (2.1)  into  (2.2)  gives: 

V24>  -  0,  (2.3) 

that  is  the  solution  of  the  irrotational  motion  of  the  incompressible  fluid 
satisfies  Laplace's  equation.  This  equation  is  to  be  solved  on  the  boundary 
z»0,  subject  to  application  of  an  axisymmetric  disturbance. 


FIGURE  2  Physical  configurations  for  a  cylindrical 
volume  of  water  (w)  moved  rigidly  upwards. 


Consider  the  case  of  a  body  of  fluid  with  air  above  and  describe  the 
interface  or  the  free  surface  by: 
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i^.v. 


f(r,z,t)  -  n(r,t)  -  z  =  0 


(2.4) 


neglecting  surface  tension,  f  is  a  function  which  is  always  zero  for  a 
particle  on  the  free  surface.  Any  disturbance  of  the  surface  clearly  implies 
some  motion  of  the  air  and  here  we  assume  that  the  pressure  in  the  air  due  to 
the  motion  is  negligible  and  the  air  pressure  may  be  approximated  by  its 
undisturbed  value.  The  boundary  conditions  at  the  undisturbed  free  surface  in 
the  context  of  the  linear  theory  can  be  written  as  (Whitham  (1974)): 


and 


on  z  =  n(r,t). 


(2.5a) 


♦t  +  gn  -  o, 


(2.5b) 


where  g  is  the  gravitational  acceleration.  Differentiation  with  respect  to  t 
of  the  kinematic  free  surface  condition  (2.5b)  and  substituting  (2.5a)  into 
this  expression,  yields  the  following  equation  to  be  satisfied  on  the  boundary 
z  =  0: 


$  +  g$  *  0  .  (2.6) 

tt  z 

On  the  rigid  boundary  z  *  -hQ  the  normal  velocity  of  the  fluid  vanishes: 

n.V<j>  =  0,  (2.7) 


where  n^is  a  vector  normal  to  the  surface.  Thus  the  problem  is  to  find  a 
solution  of  Laplace's  equation  subject  to  the  boundary  conditions  (2.5)  and 
the  initial  conditions: 


<|>(rf0,0)  -  0, 

and 

4>t(r,0,0)  =  0. 


(2.8a) 

(2.8b) 


If  hQ  is  the  depth  below  the  water  surface  at  which  the  fluid  disturbance 
originates,  then  a  further  condition  on  the  velocity  potential  is: 


r  <  r  ,  z  ■  -h  ,  (2.9) 

o  o 


where  the  function  K ( t)  is  related  to  the  generating  mechanism  of  the 
disturbance.  In  order  to  derive  closed  form  solutions  it  is  necessary  to 
choose  €(t)  appropriately,  and  this  is  further  discussed  in  $  3.0.  Once  $  has 
been  found  subject  to  the  above  conditions  the  fluid  surface  elevation  z(r,t) 
is  given  by: 


z(r,t)  -  n(r,t)  -  -  -  $  (r,o,t). 

9  t 


(2.10) 
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3.  SOLUTION  BY  METHOD  OF  INTEGRAL  TRANSFORMS 


A  solution  to  Laplace's  equation  (2.3)  subject  to  the  above  boundary 
and  initial  conditions  can  be  found  by  employing  the  standard  methods  of 
integral  transforms  (Tranter  (1971)).  In  the  case  of  the  present  problem  we 
firstly  apply  a  Hankel  transform  to  the  radial  coordinate  r,  followed  by  a 
Laplace  transform  on  the  time  coordinate  t. 

(a)  Hankel  Transform 

The  Hankel  transform  H(p)  of  a  function  H(r)  is  given  by  the 

integral: 


H(p )  =  /H(r)rJ  (pr)dr, 
n 
o 


(3.0) 


where  JR(pr)  is  the  Bessel  function  of  the  first  kind  of  order  n.  The  inverse 
Hankel  transform  is  given  by: 


H(r )  *  /H(p )pj  (pr)dp. 

n 

o 


(3.1  ) 


Laplace's  equation  (2.3)  in  cylindrical  coordinates,  the  boundary  conditions 
(2.6),  the  initial  conditions  (2.8)  and  the  velocity  potential  (2.9)  with  n~0 
in  (3.0)  can  be  expressed  as: 


and 


zz 


P2  *  -  0, 


<t>  +  qd> 

Ttt  yYz 


0, 


$(p,0,0)  »*  4>t(p,0,0)  «  0, 


5(t)  /rJ  (pr)dr,  r  <  r  ,  z  •  -  h  , 
o  o  o 

o 

0  r  >  r 


(3.2a) 

(3.2b) 

(3.2c) 

(3. 2d) 


on  applying  the  transform  (3.0).  In  taking  the  Hankel  transform  of  (2.3)  it 
has  been  assumed  that  JQ(pr)  satisfies  the  differential  equation: 

J"(pr)  +  (pr)  'j'ipr)  +  J  (pr)  *  0, 
o  o  o 

where  (  ) '  denotes  differentiation  with  respect  to  pr. 


(b)  Laplace  Transformation 


If  we  denote  by  S  (a)  the  Laplace  transform  of  a  function  S(t),  so 


that 


S  (a)  *  / S(t)e  atdt, 


(3.3) 


then  the  operation  of  multiplying  (3.2)  by  the  kernel  of  the  Laplace  transform 
and  integration  with  respect  to  t  between  o  and  •  leads  to: 


—  *  2  —  * 

*zz  “  p  ♦  "  °' 

2  —  *  —  * 

a  A  ♦  g$  ■  0;  on  *  *  0, 
z 


(3.4a) 

(3.4b) 
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—  *  .f  * 

<{>  ■  p  r  £  (o)J,  (pr  );  r<r  on  z  *  -  h  , 

Tz  o  1  o  o  o 


(3.4c) 


where  the  initial  conditions  have  been  invoiced  in  talcing  the  Laplace  transform 
of  (3.4b).  Furthermore  the  integral  property  of  Bessel  functions: 


/rVJv_1(r)dr  *  rVJv(r)  +  c, 


has  been  used  to  rewrite  (3.4c)  in  non-integral  form.  The  problem  is 
therefore  reduced  to  finding  a  solution  of  (3.4a)  subject  to  (3.4b)  and 
(3.4c). 

The  appropriate  solution  of  (3.4a)  which  satisfies  (3.4b)  and  (3.4c)  is: 


-  r  £  (a) 
o 


[a  sinhpz  -  gpcoshpz)  (pr  ) 

2  2  ™ 
p  tsinhph  +  a  coshph  ) 
o  o 


(3.5) 


and  for  waves  propagating  on  the  surface, (z  *  0)  (3.5)  can  be  expressed  as: 


4  (P,o,u) 


-gr  £  (a)J  (pr  ) 
O  I  o 


ptgsinhph  +  a  coshph  ] 
o  o 


(3.6) 


The  conversion  of  (3.6)  to  4  is  accomplished  by  applying  the  formula  (3.1): 


4  (r ,o,o)  *  -  grQ£ 


*  C  1 

(o)  /  - 

ft  Frto  S 


J, (pr  )J  (pr)  dp 
o  o 


(3.7) 


p [gsinhph  +  u  coshph  ) 
o  o 


The  inversion  of  (3.7)  to  4(r,o, t)  is  now  accomplished  by  the  use  of  the 
inverse  Laplace  transform.  Writing  (3.7)  as: 


*  *  *  , 

4  (r,o,a)  *  £  (a)B  (p|e). 


(3.8) 


where 


8*(p|a):  -  gr  f- - - 

ft  fne  i  n 


J,  (pr  T  (pr)  dp 
o  o 


(3.9) 


ptgsinhph  +  a  coshph  ] 
o  o  o 


and  applying  the  Faltung  theorem  for  Laplace  transforms,  the  velocity 
potential  is  given  by 

t 

4 (r,o, t)  -  /£( t  -  X)S(p|X)dX,  (3.10a) 

o 

or  alternatively, 

t 

4 (r,o,t)  -  /£(X)B(p|t  -  X)dX.  (3.10b) 

o 

Substitution  of  (3.10a)  (or  (3.10b))  into  (2,10)  gives: 

z(r, t)  -  -  -  3  4(r,o,t)  -  -  -  f£(X)3  B(p|t  -  X )dX  +  £{t)S(p|o),  (3.11) 

g  t  g  ^  1 

for  the  free  surface  elevation.  The  partial  derivative  Bt(p|t  "  X)  appearing 
in  (3.11)  can  easily  be  evaluated  and  is: 
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/•  cos (/pgtanhph  (t  -  X))dp 

3.B(p|t  -  X)  =  -  gr  /J,  (pr  )J  (pr) -  ° 

t  1  o  J  i  o  o 


coshph 


(3.12) 


Writing  this  integral  as  &  (p|t  -  X)  and  substitution  into  (3.11)  yields: 

z(r,t)  =  -  -  /C(X)8+(pjt  -  X)dX  +  £(t)B(p|o), 
g  o 

or  from  a  change  in  variable: 


z(r,t)  *  -  -  JC< t  -  X )8+(p j  X )dX  +  5(o)B(p|t).  (3.13) 

9  o 

Since  {  is  by  definition  the  response  of  the  generating  mechanism  of  the 
disturbance,  which  is  initially  passive,  it  follows  that  £(o)  =  0.  The  second 
term  appearing  on  the  right  side  of  (3.13)  can  therefore  be  neglected. 


In  general,  integrals  of  the  form  (3.13)  are  intractable  even  for 
simple  forms  of  £ (t  -  X)  .  If  solutions  are  to  be  found  it  is  necessary  to 
resort  to  either  asymptotic  analysis  or  numerical  integration. 


One  solution  to  (3.13)  can  be  found  if  we  suppose  that  the  fluid  is 
initially  at  rest  in  the  region  |rj<  a,  z  =  -  hQ  and  that  at  t  *  0  an 
hydrodynamic  axisymmetric  impulse  is  delivered  which  generates  the  velocity 
potential  (2.9).  If  this  impulse  acts  only  briefly  then  (2.9)  can  be  replaced 
by  a  Dirac  delta-function,  that  is  we  replace  £ (£  -  X )  by  6 (t  -  X )  in 
(3.13).  Now  (3.13)  can  be  integrated  provided  B  is  continuous  at  t  *  X  and 
o<X<t.  This  leads  to  the  result: 

z (r, t)  -  -  B+(p|t), 

9  1 

cos ( t/ pgtanhph  )dp 

- — - - —  .  (3.14) 

coshph 


or 


z(r, t) 


JJ1 


(pr  )J  (pr) 
o  o 


4.  ASYMPTOTIC  BEHAVIOUR 


The  integral  in  equation  in  (3.14)  can  be  evaluated  on  applying  the 
Kelvin  principle  of  stationary  phase,  Whitham  (1974).  Firstly  however  (3.14) 
must  be  reduced  to  canonical  form. 

fluctuating  function  which  has  an  asympototic  expansion: 

v2 


For  large  values  of  a,  JQ(a)  is  a 


v>  -  fcs> 


cos (a  -  —) , 


(4.1  ) 


so  for  large  values  of  r,  equation  (3.14)  behaves  asymptotically  as: 
“  J. (pr  ) 


ro  C  1  (pro'  r 
2(r’t}  ~  7W*  J  /pcoshph"  {e 


i (-t/pgtanhph  +  pr  -  x/4) 


-i  (t/pgtanhph  +  pr  -  ir/4).  . 

+  e  o  )  dp  , 
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and  contains  both  outgoing  waves  and  incoming  waves.  Since  r>o,  the  main 
contribution  in  the  case  t>o  comes  from  the  first  term  and  represents  waves 
propagating  outward  from  the  disturbance  to  infinity.  For  these  waves  the 
integral  can  be  expressed  as: 


z  (r,  t) 


r  r 

o  i  1 

/itr  j  /p co 


J„ (pr  )e 
o 


-iX(p  )t 


/ pcoshph 


dp, 


where 


X (p ) : =  W(p)  -  prt 


-1 


(4.2) 

(4.3) 


The  frequency  relation: 


W(p):*  / pgtanhph 

o 


(4.4) 


is  the  usual  dispersion  relation  for  gravity  waves  of  the  differential 
equation  (2.3),  where  p  is  the  wave  number  and  has  the  following  algebraic 
properties.  If  p>0,W'(p)  is  a  monotonic  decreasing  function  and  W"(p)<0.  The 
waves  propagate  to  the  right  with  the  phase  speed: 


c(p) 


W(p) 

P 


(4.5) 


and  group  velocity 


C(p)  =  D  W(p )  ••  V2  c(p)(1  +  2ph  cosech(2ph  ))  .  (4.6) 

p  00 

For  the  motions  we  are  interested  in,  that  is  the  behaviour  for  large  r  and  t, 
the  interesting  limit  is  t  +  •  with  rt-1  held  fixed.  With  (4.2)  in  canonical 
form  for  applying  Kelvin's  method  of  stationary  phase  the  surface  elevation  is 
given  by: 

r  (prt jw" (p ) | )”  ^ 

z(r,t) - ^-r - r-r - J,  (pr  )  cos  (t/ pgtanhph  -  pr),  (4.7) 

2coshph  1  o  o 

o 

where  p(r,t)  is  the  positive  root  of 


X'  (p)  =  W'(p)  -  «  0,  p>0,  >0. 

Scaling  the  independent  variables  as: 

*  -t  *  *  -1  V> 

p  ♦  k  h  ,  r  ♦  r  h  ,  t  ♦  t  (h  g  )  2 
000 


(4.8) 


(4.9) 


the  displacement  of  the  surface  can  finally  be  expressed  as, 

W'(k)  \1/2r 

z(r,t)  ~ 


1  k|w"(k)  ' 


.  J, (kr  /h  ) 
o  1  0  0 


2h  rcoshk 
o 


cos (t/ktanhk  -  kr) 


(4.10) 


after  dropping  tl^e  *.  Since  surface  water  waves  can  travel  at  a  maximum 

velocity  of  (ghQ)^2the  water  surface  remains  stationary  until  such  time  as  a 

wave  arrives  and  this  event  occurs  when  t  **  r//gh  .  If  t<r//gh  then  the 

o  o 

surface  is  stationary  and  z(r,t)  *  0,  whilst  if  t>r/rghQ,  z(r,t)  is  given  by 
(4.10). 
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Ir>  terms  of  (4.3)  and  (4.4)  equation  (4.10)  can  be  written 

z(r,t)  =  Re{A(r,  t)eiX^ ,  (4.11) 


where  the  amplitude  is: 

(W  (k) 
k|w" (k) 

The  oscillating  wave  train  represented  by  (4.11)  has  A,  k  and  W  as  functions 
of  r  and  t.  In  order  to  delineate  graphically  the  wave  train  we  take  a 
particular  value  of  kQ.  Thus 

r  =  /gh  W (k  )t, 
o  o 

represents  an  observer  moving  with  a  velocity  W'(kQ)  who  observes  a 
propagating  wave  with  wave  number  kQ  and  frequency  W(kQ).  If  ro>hQ,  the  wave 
profile  is  shown  in  Figure  3  a,b.  The  wave  profiles  are  non-periodic  with 
amplitudes  decreasing  to  zero  when  either  J1 (krQ/h0)  =  0  or 
cos (t/ktanhk  -  kr)  *  0,  that  is,  when 


V/2  r  J  (kr  /h 
\  o  1  O  O 

[" I  2hQrcoshk  * 


r  k 


where  j,  g  are  the  zeros  of  J1  (i.e.,  j1  1 


3.83171,  j12  -  7.01559...) 


kr  +  ir  (n  +  1/2) 
or  t  =  . . . 1 — 

/ktanhk 


n  **  0,1,2,  ... 


Furthermore  the  amplitude  of  the  group  decreases  as  time  advances  from  a 
maximum  Aq  to  zero,  then  increasing  to  a  new  maximum  A1 <Aq,  and  so  on.  The 
time  interval  between  wave  crests  at  earlier  times  is  longer  than  at  latter 
times  and  the  wave  profile  is  equivalent  to  the  one  given  in  section  1,  and 
represents  a  tsunami  (Adams  (1970)). 
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FIGURE  3a  A  section  of  the  wave  profile  generated  by 
a  underwater  disturbance  at  r*2000  for  r  “50,  h  “10. 
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CONCLUSIONS 


The  theoretical  results  shown  in  this  report  simulating  tsuanami 
waves,  and  the  observations  made  at  the  Bikini  nuclear  test  site  are  in  close 
agreement  within  the  context  of  the  linear  theory  of  waves.  The  primary 
assumption  employed  in  the  theoretical  analysis,  that  the  disturbance  can  be 
represented  by  a  Dirac-delta  function,  appears  to  be  acceptable.  Moreover  the 
disturbance  is  represented  mathematically  by  a  vertical  piston  type  action, 
where  the  piston  represents  the  driving  mechanism  of  the  disturbance.  This 
assumption  is  necessary  if  the  analysis  is  to  be  tractable.  Although  the 
analysis  employs  asymptotic  approximations  the  wave  structure  is 
representative  of  that  given  in  Section  1  if  rQ>h0. 

The  analysis  presented  herein  offers  some  insight  into  the  structure 
of  tsuanami  waves  generated  as  a  consequence  of  either  (i)  a  thermonuclear 
detonation,  (ii)  submarine  volcanic  explosion;  or  (iii)  abrupt  uplifting  of  a 
section  of  the  seabed.  The  magnitude  of  the  disturbance  can  be  determined 
from  the  generating  mechanism,  the  dimension  of  the  area  of  origin,  and  depth 
and  displacement  of  water. 
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